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1  Abstract 


In  this  project  we  aim  to  construct  a  high  fidelity  boundary  condition  module  for  Maxwell’s  equations 
that  can  be  interfaced  with  major  time-domain  electromagnetics  solver  systems.  There  is  ample  need 
in  the  EM  modeling  community  for  reliable  and  stable  far  field  boundary  conditions  of  high  accuracy. 
Most  existing  methods  are  limited  in  one  or  more  of  these  requirements,  and  recent  developments  in 
the  CRBC  procedure  (as  originally  presented  by  Hagstrom  and  Warburton  in  2009),  have  made  the 
technique  an  attractive  candidate  for  implementation  in  multi-purpose  solvers.  In  phase-I  of  this  project 
we  implemented  and  improved  upon  many  aspects  of  this  method,  particularly  in  light  of  the  needs  of 
high  order  accurate  Maxwell  equations  solvers  (based  on  the  discontinuous  Galerkin  method).  Error 
bounds  were  computed  and  demonstrated  for  a  number  of  cases.  We  continue  in  the  second  phase  of 
this  project  to  improve  upon  the  robustness  of  this  method,  as  we  develop  a  software  platform  which 
shall  be  its  flagship  (and  open  source)  implementation.  In  this  second  quarterly  report  we  present  a 
novel  way  to  construct  a  upwinding  numerical  flux  which  solves  the  remaining  problem  from  phase  I  of 
the  project  -  instability  issues  of  CRBC  coupling  with  DC  solver  in  2D. 

The  delay  in  submitting  this  second  quarterly  (Q2)  report  is  due  to  hire  a  new  personnel  Dr.  Ronald 
Chen  at  HyPerComp  and  having  him  come  up  to  speed  on  the  work  to  be  performed  on  this  phase 
II  contract.  From  here  on,  we  will  be  on  schedule  in  meeting  the  deliverables  (starting  with  the  third 
quartly  Q3  report  due  January  13,  2014). 
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2  Introduction 


In  this  project,  HyPerComp  is  collaborating  with  Prof.  Thomas  Hagstrom  and  his  research  group  at 
the  Southern  Methodist  University  (SMU).  Roles  of  the  two  organizations  are  very  broadly  divided 
into  mathematical  method  development  (led  by  SMU)  and  implementation,  software  development  and 
maturation  (led  by  HyPerComp).  The  project  is  coordinated  via  a  series  of  in-person  and  telephone 
meetings.  We  have  been  conducting  weekly  telephone  meetings.  Two  students,  John  Lagrone  and  Fritz 
Juhnke  have  been  included  in  the  team  and  have  been  actively  participating  in  the  work  so  far. 

Tasks:  The  following  is  a  list  of  tasks  to  be  performed  in  this  project. 

1.  Project  Formulation 

2.  Software  Development 

3.  Verification  &  Validation 

4.  Coupling 

5.  Efficiency  Testing 

6.  Release  of  software 

7.  Documentation 

8.  Sustainability  Plan 

9.  User  Support 

At  present,  we  are  working  on  solving  the  remaining  instability  problem  from  phase  I  and  testing  it 
in  sample  problems.  Primary  concerns  pertaining  to  method  stability  at  corners,  particularly  in  3D  are 
being  addressed.  CRBC  implementations  in  finite  difference  schemes,  DC  (in  FORTRAN  as  well  as  in 
MATLAB)  are  available  from  prior  research  in  this  project,  for  testing. 

We  are  presently  aiming  to  integrate  the  CRBC  module  with  the  following  codes: 

•  HDphysics  from  HyPerComp,  a  high  order  DC  based  solver 

•  MEEP  from  MIT,  an  open  source  FDTD  code 

•  cgmx  part  of  “Overture”  suite  of  simulation  codes  from  LLNL  -  high  order  finite  differences, 
second  order  PDEs 

•  CLAWPACK  a  finite  difference  suite  of  solvers  from  U. Washington 

Students  from  SMU  focus  on  an  FDTD  implementation  of  CRBC  in  the  Yee  scheme  for  a  3D 
waveguide.  Thomas  Hagstrom  starts  to  fomulate  the  CRBC  for  Maxwell’s  equation  in  3-|-I-Dimensions 
include  edges  and  corners.  We  are  in  the  process  of  developing  software  requirements  for  each  of 
the  systems  mentioned  above,  so  that  we  can  outline  a  common  implementation  of  the  method  and 
programming  techniques.  We  will  also  begin  to  implement  prototype  of  CRBC  with  DC  for  Maxwell’s 
equation  in  3D.  These  shall  be  discussed  in  the  forthcoming  report. 
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3  Instability  issue  solved  -  upwinding  flux  on  CRBC 


In  phase  I  of  the  project,  we  have  experienced  some  instability  issue  with  the  CRBC  coupling  with  DC 
when  Ni,c  =  3  and  N^g  >  10  on  unstructured  meshes  (We  have  never  encountered  any  stability  issus  on 
the  structured  mesh).  Since  we  didn’t  know  what  is  an  appropriate  upwinding  numerical  flux  to  use,  we 
simply  choosed  central  flux  for  CRBC.  This  choice  turns  out  to  be  instable  even  after  we  reformulate  the 
CRBC  on  the  corner  and  forbid  to  split  the  element  adjacent  to  the  corner  (See  next  section  for  more 
details  of  the  recap  of  the  instability  issue).  After  more  comprehensive  tests  and  debug,  we  have  found 
that  unbalanced  numerical  flux  between  volume  (upwinding  flux)  and  CRBC  boundary  (central  flux) 
caused  the  instability,  since  there  isn’t  enough  dissipation  for  the  CRBC.  Now  with  the  new  Maxwell 
like  upwinding  numerical  flux  term,  we  are  able  to  solve  the  instability  issue.  A  rigorous  proof  of  the 
stability  have  not  been  achieved  yet.  However,  we  haven’t  experienced  any  instability  issue  in  numerous 
testing  problems  with  varying  Nbc,  N^g,  dt  (See  next  section).  Details  of  mathematical  definitions  of 
the  numerical  flux  terms  are  described  below. 


Consider  the  TM  Maxwell  system: 


dE^ 

~dt 


dH^  1  dE^ 
dt  /i  dy 

dny  1  dE^ 

dt  y,  dx 

1  dny  1  dH^ 

- + - 

e  dx  e  dy 


0 

0 

0 


(1) 

(2) 

(3) 


and  set  c  =  .  Consider  a  portion  of  the  radiation  boundary  with  unit  normal  n  pointing  outward 

from  the  computational  domain  and  unit  tangent  vector  t: 


n  = 


(4) 


=  E]±yJ^  {-nyH^  +  n^Ry)  , 

From  the  CRBC  iteration  of  TM  Maxwell  equation,  we 
boundary. 


Hr,,g=n,HJ  +  nyHy,  (5) 

are  solving  a  Id  system  of  PDFs  along  the 


n  A  \  ^  1  Sin^dj  n  1  ^  rr 

(  cos  (t>g )  ^  ^ 

(I  +  cos^j)^R-j  + 

•’  dt  T  cos  4>j  e  dr 


n  A  \  ^  U  1  sin  1  5  rr 

(I  +  coscpj)  R+  j-i  ™  ,  R+  j-i  ...  TIji  j-\ 

■I'  dt  T  cos (j)j  e  dr 


T  ,  d 


{l-cos4>g)—R+y-- 


1  sin^  (j) 


1  d 


T  COB(j)j 


^  R+y  +  -^H, 


e  dr 


'n,J- 


(6) 


(7) 


— Hji  j  + - ( R+  j  +  A  —  0 

dt  2ydT 


(8) 
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This  yields  the  local  semidiscrete  scheme  equation  (9),  (10),  (11) 

1 


j  = 

dt 


(1  +  COs4>j) 
1  sin^  0 


d 

{(l-COS(/)j)— 


'o  1  sin^  (bj 

it- .j  ~  - ; —  -it- j - 1 


T  cos(/)j  ’  T  cos  (f)j 

+  \Dr{Hn,j-i-  Hn,j)+\{JMy^  (f  n  ■  (fi  -  y)l{x)dx]} 
e  A  Jx, 


Ir  1 

dt  ^  (1-t-cos^j) 


-  d 

{{l-cos(j)j)—R+j 


1  sin^  (j)j 


1  sin 


-^+,7-1 


T  cos(j)j  T  cos(j)j 

-  Hnj-i) +  1-{JM)~^  (f  n  ■  {f2  -  f2)l{x)dx]} 
e  Z  Jxi 


(9) 


(10) 


^-ff«,j  =  -^{-D^(i?+,j  +  -R-.i)  +  ^(JM)  ^  n-ifs- f3)l{x)dx}  (11) 

The  numerical  flux  for  each  iteration  are 

n-(/i-/i*)  =  n-([[-ff„.,]]-[[H„,,_i]]) 

-  -  [[^-o-i]]  +  [[^.o-]]  -  [[^.o-i]])  (12) 

n-(/2-/2)  =  n-([[if„.,_i]]-[[-ff„,,]]) 

+  -  [[!?-, ,-i]]  +  [[I?..,]]  -  [[l?.o-i]])  (13) 

n-ih-f;)  =  n-  ([[i?,.,]]  +  [[i?_,,]])  -  (14) 

[[?]]  =  9”  “  ^  =  “1)  1  depending  on  the  direction  since  we  are  solving  a  Id  system  of  PDEs.  Note 

that  we  are  solving  i?_  j  in  an  increasing  order  and  R+j  in  a  decreasing  order. 
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4  Numerical  validation 

4.1  Recap  of  Instability  issue  in  phase  I 

The  original  implementation  of  CRBC  (using  central  flux)  coupling  with  DG  on  a  unstructured  grid  is 
not  completely  stable,  especially  for  some  combination  of  DG  order  N^g,  CRBC  order  and  timestep 
size  dt  (for  example,  Ndg  =  10,  Nbc  =  3).  After  many  tests  on  different  grids  (see  figures  1,  2  (a)(b)) 
with  different  topology,  grid  size,  domain  and  time  step  size,  we  confirm  that  it  will  become  unstable 
eventually.  The  fix  in  phase  I  by  enforcing  only  one  element  along  the  corner  is  not  sufficient.  The 
only  major  difference  between  mesh  figures  1  (c)  and  (d)  are  at  the  corners,  and  in  both  cases  the 
solution  start  to  blow  out  around  time  T  =  15  with  Finaltime  =  100.  From  figure  3,  4  (little  white 
triangular  sharped  noise  along  the  boundary),  we  can  see  that  instability  starts  between  the  CRBC 
boundary  elements  and  accumulates  along  the  boundaries,  which  is  independent  of  whether  or  not  the 
corners  have  been  cutted,  nor  the  size  of  the  domain.  With  default  timestep  size  setting,  Ndg  =  10 
and  Nbc  ^  3,  in  every  case  the  solution  starts  to  blow  up  around  time  T  =  10,  Finaltime  =  100,  except 
for  grid  squarecartnx50(see  figure  1(e))  when  all  the  triangular  elements  are  obtained  by  splitting  the 
square  elements.  However,  this  still  become  unstable  when  much  smaller  timestep  is  used  (see  figure  5, 
dt  =  1.25e  -  3  is  default  timestep  size). 

The  reason  is  that  there  is  not  enough  dissipation  to  keep  the  method  stable  when  central  flux  is 
used  along  CRBC.  The  upwinding  numerical  flux  we  used  on  volume  integration  to  keep  the  method 
stable  does  not  provide  enough  dissipation  for  the  CRBC  to  be  stable.  And  there  is  no  dissipation 
coming  from  the  boundary  since  only  central  flux  is  used  on  CRBC.  The  smaller  timestep  size  is,  the 
less  dissipative  the  method  becomes.  From  figure  6,  it  is  more  obvious  to  see  how  different  timestep  size 
affect  the  dissipation  of  the  original  algorithm.  The  instability  is  not  very  strong  and  vanishes  when  Nbc 
is  big  even  if  the  central  flux  is  used  on  CRBC.  Empirically,  when  Nbc  =  3,  provide  enough  dissipation 
for  the  method  to  be  stable  with  Ndg  <  5,  and  when  Nbc  ^  4,  it  provide  enough  dissipation  for  at  least 
Ndg  =  10.  Practically  speaking,  even  if  we  use  central  flux  on  CRBC,  when  Nbc  is  big,  we  usually  don’t 
see  any  instability,  but  smaller  timestep  size  will  reduce  the  stability  of  CRBC. 
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(a)  0125 


(b)  00625 


(c)  free4 


- 

- 

-1  -0.5  0  0.5  1 

(e)  squarecartnx50 


(d)  free4goodcorner 


(f)  free4wide 


Figure  1:  List  of  the  meshs  part  I. 
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Figure  2:  List  of  the  meshs  part  II. 


4.1  Recap  of  Instability  issue  in  phase  I 


4  NUMERICAL  VALIDATION 


(a)  id=18 


(b)  id=324 


(c)  id=364 


(d)  id=424 


Figure  3:  Snapshot  of  the  instability  along  the  boundary  at  4  different  time  for  grid  free4goodcorner(see 
figure  1(d))  using  old  flux  term. 


9  of  13 


4.1  Recap  of  Instability  issue  in  phase  I 


4  NUMERICAL  VALIDATION 


Result 

10 

0.01 

0.0001 

le-6 


(a)  id=137 

Result 

10 

;1 

! 

0.01 

0.0001 

le-6 

(b)  id=207 

Result 

10 

:1 

0.01 

0.0001 

le-6 

(c)  id=247 

Result 

10 

:1 

0.01 

0.0001 

le-6 

(d)  id=287 

Figure  4:  Snapshot  of  the  instability  along  the  boundary  at  4  different  time  for  grid  longchannel(see 
figure  2(a),  2(b))  using  old  flux  term. 
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quad32  grid,  original  dt=1.1126e-3 


Time 


Figure  5:  Error  on  the  boundary  nodes  for  grid  quad32(see  figure  1(e))  with  different  timestep  size 
using  old  flux  term. 


0625  grid,  original  dt=8.2134e-4 


Time 


Figure  6:  Error  on  the  boundary  nodes  for  grid  00625(see  figure  1(b))  with  different  timestep  size  using 
old  flux  term. 
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4.2  Unconditionally  stable  with  upwinding  numerical  flux  in  2D 

The  instability  issue  is  completely  solved  in  our  new  implementation  of  the  CRBC  with  new  upwinding 
numerical  flux  along  the  CRBC  boundary  as  described  in  last  section.  We  have  tested  the  new  method 
on  all  meshes  (see  figure  1,  2)  with  various  mesh  size,  topology,  geometry,  computational  domain,  time 
step  size  dt,  DC  order  DC  Ndg,  CRBC  order  7V{,c,  and  they  are  all  stable  upto  time  finaltime  T  =  100. 
Details  of  numerical  experiments  are  summarized  in  table  4.2.  Error  results  along  time  are  shown  in  7, 
which  clearly  shows  great  advantage  over  our  old  implementation  of  CRBC  for  TM  Maxwell’s  equation 
in  2D. 


gridname 

Ndg 

Nbc 

K 

orig.  dt 

running  dt 

Note 

0125 

20 

3 

568 

5.36E-004 

l.OOE-004 

square  shape 

00625 

10 

3 

2310 

8.21E-004 

l.OOE-004 

00625 

10 

3 

2310 

8.21E-004 

5.00E-004 

square_cart_nx50 

10 

3 

5000 

7.12E-004 

5.00E-004 

freed 

10 

3 

2890 

6.93E-004 

l.OOE-004 

freed 

10 

3 

2890 

6.93E-004 

5.00E-004 

free4_goodcorner 

10 

3 

3198 

6.00E-004 

5.00E-004 

free_regoct_2 

15 

3 

2210 

3.57E-004 

2.50E-004 

octagon 

10 

3 

7.26E-004 

5.00E-004 

10 

5 

7.26E-004 

5.00E-004 

10 

9 

7.26E-004 

5.00E-004 

8 

5 

l.lOE-003 

5.00E-004 

8 

9 

l.lOE-003 

5.00E-004 

long_channel 

10 

3 

16658 

l.lOE-003 

5.00E-004 

long  channel 

10 

9 

l.lOE-003 

5.00E-004 

16 

16 

4.67E-004 

4.55E-004 

handdraw 

10 

3 

1823 

2.30E-004 

2.00E-003 

arbitrary  polygon 

free_4_wide 

10 

3 

4048 

6.27E-004 

5.00E-004 

rectangle 

Table  1:  CRBC  with  new  flux  term  experiments  summary. 
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4.2  Unconditionally  stable  with  upwinding  numerical  Bux  in  2D  4  NUMERICAL  VALIDATION 


(a)  Error  for  handdraw  polygon  and  a  rectangle. 


(b)  Error  for  32  by  2  long  channel. 


laps. 


Figure  7:  Error  plot  for  various  meshes  and  parameters  using  new  flux  term  in  CRBC. 
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